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^ ■ Abstract 

^ ' The maximum density still life problem (MDSLP) is a hard constraint optimization 

, problem based on Conway's game of life. It is a prime example of weighted constrained 

optimization problem that has been recently tackled in the constraint-programming com- 
munity. Bucket elimination (BE) is a complete technique commonly used to solve this 
kind of constraint satisfaction problem. When the memory required to apply BE is too 
high, a heuristic method based on it (denominated mini-buckets) can be used to calculate 
' bounds for the optimal solution. Nevertheless, the curse of dimensionality makes these 

, techniques unpractical for large size problems. In response to this situation, we present 

' a memetic algorithm for the MDSLP in which BE is used as a mechanism for recom- 

O . bining solutions, providing the best possible child from the parental set. Subsequently, 

a multi- level model in which this exact/metaheuristic hybrid is further hybridized with 
branch-and-bound techniques and mini-buckets is studied. Extensive experimental results 
^ ' analyze the performance of these models and multi-parent recombination. The resulting 

, algorithm consistently finds optimal patterns for up to date solved instances in less time 

than current approaches. Moreover, it is shown that this proposal provides new best known 
solutions for very large instances. 
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1. Introduction 



The game of life was proposed by John H. Conway in the 60s. Afterwards, it was divulged 
by Martin Gardner in his Scientific American columns (Gardner, 1970, 1971, 1983). The 
^ . game is played on an infinite checkerboard in which the only player places checkers on some 

H ' of its squares. Each square on the board is called a cell and has eight neighbors: the eight 

cells that share one or two corners with it. A cell is alive if there is a checker on it, and is 
dead otherwise. The contents of the board evolve iteratively, in such a way that the state 
at time t determines the state at time t + 1 according to three simple rules (see Fig. 1): 

1. If a cell has exactly two living neighbors, then its state remains the same in the next 
iteration. This is called the life constraint. 

2. If a cell has exactly three living neighbors, then it is alive in the next iteration. This 
is called the birth constraint. 

3. If a cell has fewer than two or more than three living neighbors, then it is dead in 
the next iteration. These are called the death by isolation and death by over-crowding 
constraints respectively. 
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(c) Death by isolation (d) Death by over-crowding 

Figure 1: Rules for the Game of Life 

As it can be seen, the game of hfe is defined in terms of simple rules, but these can never- 
theless generate incredibly complicated patterns and dynamics, and hence, it has attracted 
the interest of many scientists. 

One challenging constraint optimization problem based on the game of life is the maxi- 
mum density still life problem (MDSLP). In order to introduce this problem, let us define a 
stable pattern (also called a still life) as a board configuration that does not change through 
time, and let the density of a region be its percentage of living cells. The MDSLP in an n x n 
grid consists of finding a still life of maximum density. Elkies (1998) has shown that, for 
infinite boards, the maximum density is 1/2 (for finite size, no exact formula is known). In 
this paper, we are concerned with the MDSLP and finite patterns, that is, finding maximal 
n X n still lifes. 

The MDSLP is very hard to solve, and though it has not been proven to be NP-hard to 
the best of our knowledge, no polynomial-time algorithm for it is known. Our interest in 
this problem is manifold. Firstly, it must be noted that the patterns resulting in the game 
of life are very interesting. For example, by clever placement of the checkers and adequate 
interpretation of the patterns, it is possible to create a Turing-equivalent computing machine 
(Berlekamp, Conway, & Guy, 1982). From a more applied point of view, it is interesting to 
consider that many aspects of discrete dynamical systems have been developed or illustrated 
by examples in the game of life (Gardner, 1971, 1983). In this sense, finding stable patterns 
can be regarded as a mathematical abstraction of a standard issue in discrete systems 
control. Finally, the MDSLP is a prime example of weighted constrained optimization 
problem. As such, it constitutes an excellent test bed for different optimization techniques. 
Indeed, the problem has been included in the CSPLib^ repository. A dedicated web page^ 
maintains up-to-date results. 

The MDSLP has been tackled using different approaches. Bosch and Trick (2002) com- 
pared different formulations for the MDSLP using integer programming (IP) and constraint 
programming (CP). Their best results were obtained with a hybrid algorithm mixing the 
two approaches. They were able to solve the cases for n = 14 and n = 15 in about 6 
and 8 days of CPU time respectively. Smith (2002) used a pure constraint programming 
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Table 1: Best experimental results reported in (Bosch & Trick, 2002) (CP/IP), (Larrosa & 
Morancho, 2003) (BE) and (Larrosa et al., 2005) (HYB-BE) for solving the MDSLP. 
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approach to undertake the problem and proposed a formulation of the problem as a con- 
straint satisfaction problem with 0-1 variables and non-binary constraints. Its dual graph 
translation into a binary constraint satisfaction problem was also considered. Surprisingly, 
it was proven that, although the dual representation has as many variables as the original 
one and very large domains, its performance was much better. However, only instances up 
to n = 10 could be solved. The best results for this problem were reported by Larrosa and 
Morancho (2003) and Larrosa, Morancho, and Niso (2005), showing the usefulness of bucket 
elimination (BE), an exact technique based on variable elimination and commonly used 
for solving constraint satisfaction problems described in detail in Section 2.4. Their basic 
approach could solve the problem for n = 14 in about 10^ seconds. Further improvements 
pushed the solvability boundary forward to n = 20 in about twice as much time. Recently, 
Cheng and Yap (2005, 2006) have tackled the problem via the use of ad-hoc global case 
constraints, but their results are comparable to IP/CP hybrids, and thus lie far from the 
ones obtained previously by Larrosa et al. 

Table 1 resumes experimental results for current approaches used to tackle the MDSLP. 
The first column contains the problem size. The second column shows the optimal solution 
as the number of dead cells. Remaining columns report times in seconds by the hybrid 
IP/CP algorithm of Bosch and Trick (2002), by the BE approach of Larrosa and Morancho 
(2003) and by the BE/search hybrid of Larrosa et al. (2005). Although different computa- 
tional platforms may have been used for these experiments, the trends are very clear and 
give a pristine indication of the potential of the different approaches. In this sense, note that 
all of these techniques applied to the MDSLP are exact approaches that are inherently lim- 
ited for increasing problem sizes, and whose capabilities as anytime algorithms are unclear. 
To tackle this problem, we recently proposed the use of hybrid methods combining exact 
and metaheuristic approaches. Particularly, in (Gallardo, Cotta, &: Fernandez, 2006a) we 
considered the hybridization of BE with evolutionary algorithms (an stochastic population- 
based search method) endowed with tabu search (a local search method). The resulting 
algorithm was a memetic algorithm (MA; see Section 2.2) that used BE as a mechanism for 
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recombining solutions, providing the best possible child from the parental set. Experimental 
tests indicated that the algorithm provided optimal or near-optimal results at an accept- 
able computational cost. Afterwards, in (Gallardo, Cotta, & Fernandez, 2006b) we studied 
expanded multi-level models in which our previous hybrid algorithm was further hybridized 
with a branch-and-bound derivative, namely beam search (BS). Studies about the influence 
that variable clustering and multi-parent recombination have on the performance of the 
algorithm were also conducted. Results indicated that variable clustering was detrimental 
for this problem but also that multi-parent recombination improves the performance of the 
algorithm. To the best of our knowledge, these are the only heuristic approaches that have 
been applied to this problem to date. 

This paper includes and extends our previous research on this problem. As new contri- 
butions, we have redone all experiments using an improved implementation of the bucket 
elimination crossover operator, described in Section 3.2. Additionally, we present a more 
extensive experimental analysis of our BS/MA hybrid described in (Gallardo et al., 2006b), 
analyzing the sensitivity of its parameters. We also propose a new hybrid algorithm that uses 
the technique of mini-buckets (MB) (Dechter, 1997) to further improve the lower bounds 
of the partial solutions that are considered in the BS part of the hybrid algorithm. This 
new algorithm is obtained from the hybridization, at different levels, of complete solving 
techniques (BE), incomplete deterministic methods (BS and MB) and stochastic algorithms 
(MAs). An experimental analysis shows that this new proposal consistently finds optimal 
solutions for MDSLP instances up to n = 20 in considerably less time than all the previ- 
ous approaches reported in the literature. Finally, in order to test the scalability of our 
approach, this novel hybrid algorithm has been run on very large instances of the MDSLP 
for which optimal solution are currently unknown. Results were very successful, as the 
algorithm performed at the state-of-the-art, providing solutions that are equal or better to 
the best ones reported to date in the literature. 

The paper is self-contained and structured as follows: Section 2 gives preliminary con- 
cepts that will be used in the rest of paper. Section 3 defines the MDSLP as a weighted 
constraint satisfaction problem, and shows how to solve it using BE. In Section 4, a MA 
for the MDSLP that uses BE as a recombination operator is presented and experimentally 
evaluated along with a hybrid multilevel algorithm that integrates the previous MA with 
Branch-and-Bound derivatives. Section 5 proposes and evaluates a novel hybrid algorithm 
that exploits the technique of mini-buckets. Finally, Section 6 presents conclusions and 
outlines future work. 

2. Preliminaries 

In this section, we briefly introduce concepts and techniques that will be used in the rest of 
paper. To this end, we first present beam search, a heuristic tree search algorithm derived 
from the branch and bound method. Subsequently, memetic algorithms are introduced. 
Finally, weighted constraint satisfaction problems are defined and the technique of bucket 
elimination -commonly used to solve them- is introduced. For the sake of notational sim- 
plicity, we use in this last subsection the notation of (Larrosa & Morancho, 2003; Larrosa 
et al, 2005). 
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2.1 Beam Search 

Branch and bound (BB) (Lawler & Wood, 1966) is a general tree search method to solve 
combinatorial optimization problems. Tree search methods are constructive, in the sense 
that they work on partial solutions. In this way, tree search methods start with an empty 
solution that is incrementally extended by adding components to it. The way that partial 
solutions can be extended depends on the constraints imposed by the problem being solved. 
The solution construction mechanism maps the search space to a tree structure, in such a 
way that a path from the root of the tree to a leaf node corresponds to the construction 
of a solution. In order to efficiently explore this search tree, BB algorithms maintain an 
upper bound and estimate lower bounds for partially constructed solutions. Assuming a 
minimization problem, the upper bound corresponds to the cost of the best solution found 
so far. During the search process, a lower bound is computed for any partial solution 
generated, estimating the cost of the best solution that can be constructed by extending 
it. If this lower bound is greater than the current upper bound, solutions constructed by 
extending it will not lead to an improvement, and thus all nodes descending from it can 
be pruned from the search tree. Clearly, the capability of the algorithm for pruning the 
search tree depends on the existence of an accurate lower bound, that should additionally 
be computationally inexpensive in order to be practical. 

Beam search (BS) (Barr &: Feigenbaum, 1981) algorithms are incomplete derivatives 
of BB algorithms, and are thus heuristic methods. Essentially, BS works by extending 
every partial solution from a set B (called the beam) in at most k^xt possible ways. Each 
new partial solution so generated is stored in a set B'. When all solutions in B have been 
processed, the algorithm constructs a new beam by selecting the best up to kf,^ (called the 
beam width) solutions from B'. Clearly, a way of estimating the quality of partial solutions, 
such lower bound, is needed for this. 

One interesting peculiarity of BS is that it works by extending in parallel a set of 
different partial solutions in several possible ways. For this reason, BS is a particularly 
suitable tree search method to be used in a hybrid collaborative framework, as it can be 
used to provide periodically promising partial solutions to a population based search method 
such as a MA. Gallardo, Cotta, and Fernandez (2007) have shown that this kind of hybrid 
algorithms can provide excellent results for some combinatorial optimization problems. We 
will subsequently present a hybrid tree search/memetic algorithm for the MDSLP based on 
this idea. 

2.2 Memetic Algorithms 

Evolutionary algorithms (EAs) are population-based metaheuristic optimization methods 
inspired by biological evolution (Back, 1996; Back, Fogel, & Michalewicz, 1997). In order 
to explore the search space, the EA maintains a set of solutions known as the population of 
individuals {fi is used to denote the total number of individuals in the population). These 
are initialized usually in a random way across the search space, although an heuristic may 
also be used. After the initialization, three different phases are iteratively performed until 
a termination condition is reached: selection, reproduction and replacement. In the context 
of EAs, the objective function assigning values to each solution is termed a fitness function, 
and is used to guide the search by comparing the goodness of different individuals. 
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Memetic Algorithm 



1 : for i := 1 to popsize do 

2: pop[i] := RANDOM BOARD(n) 

3: pop[i] := LOCAL SEARCii{pop[i]) 

4: Evaluate (pop [i]) 

5 : end for 

6 : while allowed runtime not exceeded do 

7: for z := 1 to offsize do 

8: if recombination is performed {px) then 

9: parenti := SELECT(pop) 

10: parent2 := SELECT(pop) 

11: par entarity '■= Select (pop) 

12: offspring[i] := RECOMBlNE(parenti, parent2, ■ ■ ■ , parentarity) 

13 : else 

14: offspring[i] := SELECT(pop) 

15: end if 

16: if mutation is performed {pm) then 

17: offspring[i] := MlJTATE{offspring[i]) 

18: end if 

19: offspring[i] := LOCAL SEARC}i{offspring[i]) 

20: EVAhV ATE{offspring[i]) 

21 : end for 

22: pop := REPLACE(pop, offspring) 

23 : end while 



Figure 2: Pseudo code of a memetic algorithm (MA). Although different variants are possible 
with respect to this scheme, it broadly captures the algorithmic structure typically used in 
MAs. 

Note that EAs are black box optimization procedures in the sense that no knowledge 
of the problem (apart from the fitness function) is used. The need to exploit problem- 
knowledge has been repeatedly shown in theory (Wolpert & Macready, 1997) and in prac- 
tice (Davis, 1991) though (see also Culberson, 1998). Different attempts have been made 
to answer this need; Memetic algorithms (Moscato, 1999; Moscato &: Cotta, 2003; Moscato, 
Mendes, & Cotta, 2004; Krasnogor & Smith, 2005) (MAs) are one of the most successful 
to date (Hart, Krasnogor, & Smith, 2005). As EAs, MAs are also population based meta- 
heuristics. The main difference is that the components of the population (sometimes termed 
agents in the MAs terminology) are not passive entities. These agents are active entities 
that cooperate and compete in order to find improved solutions. 

There are many possible ways to implement MAs. The most common implementation 
consists of combining an EA with a procedure to perform local search on some or all solutions 
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in the population during the main generation loop (cf. Krasnogor & Smith, 2005). Fig. 2 
shows the general outline of such a MA. It must be noted however that the MA paradigm 
does not simply reduce itself to this particular scheme and there are diverse places (e.g., 
population initialization, genotype to phenotype mapping, evolutionary operators, etc.) 
where problem specific knowledge can be incorporated. In this work, apart form using tabu 
search (Glover, 1989, 1990) (TS) as a local search procedure within the MA, we have designed 
an "intelligent" recombination operator that uses an exact technique (bucket elimination) 
in order to find the best solution that can be constructed from a set of parents without 
introducing implicit mutation (i.e., exogenous information). 

2.3 Weighted Constraint Satisfaction Problems 

A weighted constraint satisfaction problem (WCSP) (Schiex, Fargier, & Verfaillie, 1995; 
Bistarelli, Montanari, & Rossi, 1997) is a constraint satisfaction problem (CSP) in which 
preferences among solutions can be expressed. Formally, a WCSP can be defined by a tuple 
{X, V, T), where V = {Di, • • • , Dn} is a set of finite domains, X = {xi, • • • , is a set of 
variables taking values from their finite domains {Di is the domain of variable Xi) and J- is 
a set of cost functions (also called soft constraints or weighted constraints) used to declare 
preferences among possible solutions. Permitted assignments of variables receive finite costs 
that express their degree of preference (the lower the value the better the preference) and 
forbidden assignments receive cost cx3. Note that each / € is defined over a subset of 
variables, var{f) C X, called its scope. The objective function F is defined as the sum of 
all functions in T: 



The assignment of value Vi € Di to variable Xi is noted Xi = Vi. A partial assignment of 
m variables is a tuple t = (xj^ = vijXi^ = V2, - ■ ■ ,Xi^ = Vm)- A complete assignment of 
all variables with values in their domains that satisfies every soft constraint (i.e., with a 
finite valuation for F) represents a solution to the WCSP. The optimization goal is to find 
a solution that minimizes this objective function. 

A WCSP'^ instance is usually depicted by means of its constraint graph, which has one 
node for each variable Xi £ X, and one edge connecting any two nodes whose variables 
appear in the same scope of some cost function f £ T . 

2.4 Bucket Elimination 

Bucket elimination (BE) (Dechter, 1999) is a generic technique suitable for many automated 
reasoning and optimization problems and, in particular, for WCSP solving. The functioning 
of BE is based upon the following two operators over functions (Larrosa et al., 2005): 



3. Observe that the constraints are not weighted in the sense of having an external weight parameter 
assigned to them. Indeed, each of them has the same influence on the overall function value, as shown in 
Eq. (f). The reason they are called "weighted" is that the output of each function is not binary (satisfied 
vs. unsatisfied) but a numerical value when it is satisfied. 




(1) 
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• the sum of two functions f and g, denoted {f + g), is a new function with scope 
var{f) U var[g) which returns for each tuple the sum of costs of / and g, 

{f + 9){t) = f{t) + g{t). (2) 

• The ehmination of variable xi from /, denoted / JJ- Xj, is a new function with scope 
var{f) — {xi} which returns for each tuple t the minimum cost extension of t to Xj, 

(/ -II Xi){t) = min^(zDi{f{t ■ {xi = v))}, (3) 

where t ■ {xi = v) means the extension of the assignment t with the assignment of 
value V to variable Xj. Observe that when / is a unary function (i.e., it has arity one), 
a constant is obtained upon elimination of the only variable in its scope. 



Bucket Elimination for a WCSP {X, V, F) 




function BE(;f, V, T) 


1 


for i := n downto 1 do 


2 


5. := {/ G ^ 1 Xi G varif)} 


3 


9i ■■= iT.feBj)^Xi 


4 


-^:= {J'[J{9^}) - B^ 


5 


end for 


6 


t :=0 


7 


for i := 1 to n do 


8 


V := argminaeDAiY^feB, /)(* " i^i = «))} 


9 


t := t ■ (xi = v) 


10 


end for 


11 


return(jF, t) 




end function 



Figure 3: The general template, adapted from Larrosa and Morancho (2003), of bucket 
elimination for a WCSP {X,V,F). 

Without loss of generality, let us assume a lexicographical ordering for the variables in 
X, i.e., o = (xi,X2,--- ,x„). Fig. 3 shows a pseudo-code of the BE algorithm for solving 
a WCSP instance, that returns the optimal cost in and one optimal assignment in t. 
Observe that, in a first phase, BE eliminates one variable Xi & X in each iteration of the 
loop comprising lines 1-5. This is done by computing firstly the bucket Bi of variable Xj as 
the set of all cost functions in having Xj in their scope. Then, a new function gi is defined 
as the sum of all these functions in Bi in which variable Xj has been eliminated. Finally, 
is updated by removing the functions involving Xj (i.e., those in Bi) and adding the new 
function that does not contain Xj. The consequence is that Xj does not exist in but 
the value of the optimal cost is preserved. The elimination of xi produces an empty scope 
function (i.e., a constant) which is the optimal cost of the problem. Then, in lines 6-10, BE 
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generates an optimal assignment of variables by considering these in the order imposed by 
o: this is done by starting from an empty assignment t and assigning to Xi the best value 
regarding the extension of t with respect to the sum of functions in Bi {argmina{f{a)} 
denotes the value of a producing minimum /(a)). 

Note that BE has exponential space complexity because, in general, the result of sum- 
ming functions or eliminating variables cannot be expressed intensionally by algebraic ex- 
pressions and, as a consequence, intermediate results have to be collected extensionally in 
tables. To be precise, the complexity of BE depends on the problem structure (as captured 
by its constraint graph G) and the ordering o. According to Larrosa and Morancho (2003), 
the complexity of BE along ordering o is time Q{Q x nx space Q{nx (°^), 

where d is the largest domain size, Q is the cost of evaluating cost functions (usually as- 
sumed 0(1)), and w*{o) is the induced width of the graph along ordering o, which describes 
the largest clique created in the graph by bucket elimination, and which corresponds to 
the largest scope of a function recorded by the algorithm. Although finding the optimal 
ordering is NP-hard (Arnborg, 1985), heuristics and approximation algorithms have been 
developed for this task (check Dechter, 1999 for details). 

3. The Maximum Density Still Life Problem 

According to the definition of MDSLP and the three rules of the game, it is easy to see that 
each cell in a still life must satisfy the following conditions: 

• If the cell is alive, it must have two or three neighbors. 

• If the cell is dead, it will have either more than three or less than three neighbors. 

Note that finite still lifes are not allowed to produce new living cells outside the grid, and 
hence stability conditions must hold in the cells surrounding the n x n square, that are 
assumed to be dead. This can equally be achieved by requiring further that: 

• if the cell is at the boundary of the nx n square, it must not be part of a sequence of 
three consecutive living cells in the direction of the boundary. 

Fig. 4 shows some maximum density still lifes for small values of n. 




Figure 4: Maximum density still lifes for n € {3,4, 5}. 

The constraints and objectives of the MDSLP are formalized in the following subsections 
in which we follow a similar notation to the one used in (Larrosa &: Morancho, 2003; Larrosa 
et al, 2005). 
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3.1 Problem Formulation 

To state the problem formally, let r be an n x n binary matrix, such that rij G {0, 1}, 1 ^ 
hj^n- {I'ij = if cell (i, j) is dead, and 1 otherwise). In addition, let Af{r,i,j) be the set 
comprising the neighborhood of cell rjj: 

A/'(r,i,j) ={r(,+^)(j-+j,)|x,y € {-1,0,1} Ax^+yVOA 1 {i + x),{j+y) ^ |Ir|| }(4) 

where ||r|| denotes the number of rows (or columns) of square matrix r, and let the number 
of living neighbors for cell r^j be noted as r]{r,i,j): 

According to the rules of the game, let us also define the following predicate that checks 
whether cell r^j is stable: 

„• .^ _ [2 < Vir,i,j) ^ 3, rij = 1 
^^"'^'^^-\ry(r,z,j)/3, r,, = 0. 

In order to check boundary conditions, we will further denote by r the (n + 2) x (n + 2) 
matrix obtained by embedding r in a frame of dead cells: 

0, otherwise. 

The maximum density still life problem for an n x n board, MDSLP(n), can now be 
stated as finding an n x n binary matrix r, such that 

(1 — rij) is minimal, (8) 



subject to 



A (9) 



3.2 The MDSLP as a Weighted Constraint Satisfaction Problem 

As shown by Larrosa and Morancho (2003) and Larrosa et al. (2005), the MDSLP fits 
nicely within the framework of WCSPs. To this end, an n x n board configuration can be 
represented by an n-dimensional vector (ri, r2, . . . , r„). Each vector component encodes (as 
a binary string) a row, so that the j-th bit of row (noted rij) signifies the state of the 
j-th. cell of the i-th row (a value of 1 represents a live cell and a value of a dead cell). 

Two functions over rows will be useful to describe the constraints that must be satisfied 
by a valid configuration. The first one, 

zeroes{a) = {1 — ai), (10) 

1 <j<n 
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returns the number of dead cells in a row (i.e., the number of zeroes in binary string a). 
The second one, 

Adjs{a) = Adjs'{a,l,0) (11) 

{I, i > n 

Adjs'{a,i + 1 ,1 + 1), CLi = 1 

max{l, Adjs'{a,i + 1 , 0)), Oj = 0, 

computes the maximum number of adjacent living cells in row a. We also introduce a ternary 
predicate, Stable{ri-i , r, rj+;), that takes three consecutive rows in a board configuration 
and is satisfied if, and only if, all cells in the central row are stable (i.e., all cells in r will 
remain unchanged in the next iteration): 

Stable{a, b, c) = Ai<.t<n '^'(a, b, c, i) (12) 



5(a, 6, c, i) 



2 ^ r/(a, 6, c, i) ^ 3, bi = 1 
^r]{a,b,c,i) / 3, bi = 

7]{a, b, C, i) = Emax(l,^-l)<:,■<min(n,^+l)(ai + + Cj) - bi, 



where 77(0, b, c, i) is the number of living neighbors of cell bi, assuming a and c are the rows 
above and below row b. 

The MDSLP can now be formulated as a WCSP using n cost functions /», i € {1 . . n}. 
Accordingly, /„ is binary with scope the last two rows of the board {var{fn) = {rn-i,rn}) 
and is defined as: 

{00, -^Stable{a, b, 0) 

00, Adjs{b) > 2 (13) 

zeroes{b), otherwise. 

The first line checks that all cells in row r„ are stable, whereas the second one checks that 
no new cells are produced below the nxn board. Note that any pair of rows representing an 
unstable configuration is assigned a cost of 00, whereas a stable one is assigned its number 
of dead cells (to be minimized). 

For i G {2 . .71 — 1}, corresponding fi cost functions are ternary with scope var{fi) = 
{ri_i, Tj, rj_|_i} and are defined as: 

{00, ^Stable{a, b, c) 

00, ai = bi = ci = l , , 

a„ = 6„ = c„ = 1 (1^) 
zeroes (b), otherwise. 

In this case, boundary conditions are checked to the left and right of the board. As regards 
cost function /i, it is binary with scope the first two rows of the board {var{fi) = {ri, r2}) 
and is specified similarly to /„: 

' CO, ^Stable{0, b, c) 

00, Adjs{b) > 2 (15) 

zeroes (b), otherwise . 



11 



Gallardo, Gotta, & Fernandez 



3.3 Solving the MDSLP with BE 

According to the formulation of the MDSLP as a WSCP introduced in Section 3.2, the 
corresponding constraint graph has a sequential structure, in which an arbitrary row is 
linked to the two rows above and below it. Due to this simple structure, it is easy to find 
an optimal elimination order for BE, and variables can be eliminated starting with the last 
one and proceeding in decreasing order. Fig. 5 shows the resulting algorithm. Function 
BE takes two parameters: n, the size of the instance to be solved, and D, the domain for 
each variable (row) in the solution. If domain D is set to {0 . . 2" — !} (i.e., a set containing all 
possible rows) the function implements an exact method that returns the optimal solution 
for the problem instance (as the number of dead cells) and a vector corresponding to the 
rows of that solution. The algorithm starts by eliminating the last variable r„, whose bucket 
is Bn = {fn, fn-i}, the Only cost functions containing r„ in their scopes. In lines 1-3, Bn is 
used to compute a new cost function, gn{a,b), with scope {r„_2,r„-i}, that represents the 
cost of the best extension of (r„_2 = a, r„„i = b) to the removed variable r„. At this point, 
the bucket of the next variable, r„_i, is Bn~i = {gn, fn~2}, that can be used to compute 
a new cost function, gn-i{a,b) with scope {r„_3,r„_2} representing the cost of the best 
extension of (r^-s = a,r„_2 = b) to the removed variables r„ and r^-i- This process can 
be iterated (lines 4-8) to eliminate variables up to r^. Optimal values for variables ri and 
r2 can be calculated using an exhaustive search (line 9). At this time, the optimal cost can 
be calculated and the optimal values for remaining variables can be set in increasing order 
using their bucket and variables assigned beforehand (lines 11-14). 

Note that the space complexity of the algorithm, when used as an exact method, is 
0(n X 2^"), due to the memory required to store extensionally n cost functions gi, having 
each 2" x 2" entries. The time complexity is 0(n^ x 2^") due to lines 4-8, as finding the 
minimum of 2" alternatives, being the computation of each one @{n), has to be repeated 
0(n X 2^") times. On the other hand, a basic search-based solution to the problem could 
be implemented with worst case time complexity 0(2'^"' •*) and polynomial space. Observe 
that the time complexity of BE is therefore an exponential improvement over basic search 
algorithms, although its high space complexity makes the approach impractical for large 
instances. 

One interesting optimization, presented by Larrosa and Morancho (2003), allows reduc- 
ing the complexity of the algorithm. In the following, we assume that n is even, although a 
similar reasoning can be used if the size of the board is odd. The optimization avoids the 
computation needed to eliminate variables ri,r2, ■ ■ ■ ,r^_i, as a result of the symmetry of 
the problem. In this way, the algorithm starts by eliminating variables r„,r„_i, . . .r^_^_2■ 
Observe that, at this point, cost functions gn,gn-i, ■ ■ ■ )5|+2 have been computed. At this 
point, the order to eliminate remaining variables can be changed to ri,r2, . . . ,r^_i. The 
elimination of ri would produce gi with scope {ri, r2}, but this computation can be avoided, 
as it is the same to eliminate ri or to rotate the board by 180 degrees and eliminate variable 



where r denotes the reflection value of the binary string r. Moreover, if the board is 
vertically reflected, an equivalent problem is obtained, so it follows that 



r. 



I.e.: 




(16) 




(17) 
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Bucket Elimination for the MDSLP 




function BE(n, V) 


1 


for a, 6 G D do 


2 


Qnia.b) := miricipvi fn-i(ci',b, c) + fn(b,c)} 


3 


end for 


4 


for i := n — 1 downto 3 do 


5 


for a,b (z do 


6 


(7, (a, b) := minc.(zT>{ fi^\ (a, b, c) + (b, c)| 


7 


end for 


8 


end for 


9 


{n,r2) := argmina,bel5{g3(a,b) +fi(a,b)} 


10 


opt := g3{ri,r2) + /i(r-i,r2) 


11 


for i := 3 to n — 1 do 


12 


ri := argmincgp{fi_i(ri„2,ri-i,c) + gi+i(ri_i, c)} 


13 


end for 


14 


r„ := argmincgp{fn_i(rn-2,rn-i,c) + fn(rn_i,c)} 


15 


return {opt, (ri,r2, . . • ,r„)) 




end function 


Figure 5: Bucket elimination for the MDSLP. 



and hence 

gi{a,b) = gn{b,a). (18) 

The optimized algorithm is obtained by applying the same reasoning to the rest of the 
variables. Larrosa and Morancho (2003) and Larrosa et al. (2005) have used this method^ 
to solve the MDSLP up to size 14. The fourth column in Table 1 reproduces their results, 
obtained with a 2GHz Pentium IV machine with 2Gb of memory. Notice the limitations of 
the approach: the n = 15 instance could not be solved due to space restraints. In Section 
4.1, we will show how BE can be embedded in a MA with reduced complexity in order to 
implement a smart recombination operator. 

4. A Multi-Level Memetic/Exact Hybrid Algorithm for the MDSLP 

WCSPs are very amenable for being undertaken with evolutionary metaheuristics. Obvi- 
ously, the quality of the results will greatly depend on how well knowledge of the problem is 
incorporated into the search mechanism. Our final goal is to present an algorithmic model 
based on the hybridization of MAs with exact techniques at two levels: within the MA 
(as an embedded operator), and outside it (in a cooperative model). Firstly, we will focus 
in the next subsection on the first level of hybridization, that incorporates an exact tech- 

4. Actually, an instance of the algorithm in which each variable is allowed to take values in the whole 
computation domain. This is not our case as it will be shown in next sections. 
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nique (namely BE) within the MA as an embedded recombination operator. Subsequently, 
we will proceed to a second level of hybridization, in which the MA cooperates with a 
branch-and-bound based beam search algorithm. 



4.1 A Memetic Algorithm with BE for the MDSLP 

In this subsection we describe a MA for the MDSLP that uses tabu search (TS) as a local 
search operator and BE as an optimal recombination operator. Before detailing these two 
components, let us describe the basic underlying evolutionary algorithm (EA). 



4.1.1 Representation and Fitness Calculation 

The natural representation of MDSLP solutions is their binary encoding. Accordingly, a 
configuration for an n x n board will be represented as a binary n x n matrix r. Clearly, 
infeasible solutions can be represented, since not all such binary matrices will correspond to 
stable patterns. One way to deal with this scenario is using penalty-based fitness functions. 
To be precise, the fitness (to be minimized) of a configuration r is defined as: 

n n n+2 n+2 

/(O = E - ^*^) + ^E E [^^M^^r^ hj)) + (1 - nj)Uv{r, hj))] ■ (19) 

i=l j=l i=l j=l 



Recall that stability is not only required within the n x n board, but also in its immediate 
neighborhood, and this is taken into account by working with r, the {n + 2) x {n + 2) binary 
matrix obtained by embedding r in a frame of dead cells, as defined in (7). K is a constant, 
r]{r,i,j) is the number of live neighbors of cell r^j, and (po, (pi : N — > N are two functions 
(to be used with dead or alive cells respectively), that take the number of alive neighbors 
of a cell, and return a penalty depending on how many of them should be fiipped to have 
a stable configuration, defined as: 




2 ^ ?7 s: 3 

i{t]) = {K' + 2-T], 77 < 2 (21) 
-3, 77 > 3, 



where K' is another constant. The first double sum in (19) corresponds to the basic quality 
measure for feasible solutions, i.e., its number of dead cells. With respect to the last term, 
it represents the penalty for infeasible solutions. The strength of penalization is controlled 
by constants K and K'. The values we have chosen for them {K = n^ and K' = bn"^) 
ensure that given any two solutions, the one that violates less constraints is preferred; 
if two solutions violate the same number of constraints, the one whose overall degree of 
violation (i.e., distance to feasibility) is lower is preferred. Finally, if the two solutions are 
feasible, the penalty term is null and the solution with the higher number of live cells is 
better. 
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4.1.2 A Local Improvement Strategy Based on Tabu Search 



The fitness function defined above provides a stratified notion of gradient that can be ex- 
ploited by a local search strategy. Moreover, notice that the function is quite decomposable, 
since interactions among variables are limited to adjacent cells in the board. Thus, when- 
ever a configuration is modified, the new fitness can be computed only considering the 
cells located in adjacent positions to changed cells. To be precise, assume that cell is 
modified in solution r, resulting in solution s; the new fitness /(s) can be computed as: 



f{s)=f{r)+K 
and functions A/i and A/2 are defined as: 



A f liuj, r]{r, i, j)) + ^h{c,r]{c) 



(22) 



'0, 



r/ = 2 



A/i(c,?7) = <^ (-1) </.o(r?), r/ = 3 



-1) 



otherwise 



A/2(c', r/, c) = (1 - c') A/2,o(r?, c) + c'A/2,i(77, c) 
A/2,o(??,c) 



A/2,1 (^,c) 



(K' + l, 


(7j = 2 A c = 0) V (r/ = 4 A c = 


1) 


-(^' + 1), 


rj = 3 




lo, 


otherwise 




(K' + l, 


(jj = 2 A c = 1) V (?7 = 3 A c = 


0) 


-(K' + l), 


(7^ = 1 A c = 0) V (r? = 4 A c = 


1) 


< 1, 


(7^=lAc=l)V(?7^4Ac = 


0) 


-1, 


(t^ = 0) V (77 ^ 5 A c = 1) 




lo, 


otherwise. 





(23) 
(24) 
(25) 

(26) 



Using this efficient fitness re-computation mechanism, our local search strategy explores 
the neighborhood N{r) = {s \ Hamming(r, s) = 1}, i.e., the set of solutions obtained by 
flipping exactly one cell in the configuration. This neighborhood comprises configura- 
tions, and it is fully explored in order to select the best neighbor. In order to escape from 
local optima, a tabu-search scheme is used: up-hill moves are allowed, and after flipping a 
cell, it is put in the tabu list for a number of iterations (randomly drawn from [n/2,3n/2] 
to hinder cycling in the search). Thus, it cannot be modified in the subsequent iterations 
unless the aspiration criterion is fulfilled. In this case, the aspiration criterion is improv- 
ing the best solution found in that run of the local search strategy. The whole process is 
repeated until a maximum number of iterations is reached, and the best solution found is 
returned. 



4.1.3 Optimal recombination with BE 

Recall that the fitness function that we have defined is able to evaluate any representable 
configuration (feasible or not), and hence, the binary representation used turns out to be 
freely manipulable. With this setting, any standard recombination operator for binary 
strings could be used in principle. For example, the two-dimensional version of single-point 
crossover (2D-SPX), depicted in Fig. 6, could be employed. Although such a blind operator 
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is feasible from a computational point of view, it would perform poorly, as it would behave 
like a macromutation operation. In order to achieve a sensible recombination of information, 
we can resort to BE. 

Even though the performance of BE as an exact method for the MDSLP was better 
than basic search-based approaches, it was shown in Section 3.3 that the corresponding 
time and space complexity were still very high, making it unsuitable for large instances. In 
the following, we explain how BE can be used to implement an intelligent recombination 
operator for the MDSLP. Such operator will explore the dynastic potential (Radcliffe, 1994) 
(possible children) of the solutions being recombined, providing the best solution that can be 
constructed without introducing implicit mutation, i.e., exogenous information (cf. Cotta & 
Troya, 2003). Moreover, we will show that this operator is tractable from a computational 
point of view. 

For this purpose, let x = (xi, a;2, • • • , Xn) and y = (yi, y2,- " > Un) be two board configu- 
rations for an nxn instance of the MDSLP. Our operator will calculate the best configuration 
that can be obtained by combining rows in x and y without introducing information not 
present in any of the parents. This can be achieved by restricting the domain of variables 
in BE to take values corresponding to the rows of the configurations being recombined. 
Using the optimized version of the BE algorithm, the recombination operator becomes 
BE-Opt(n, {xi,X2, • • • , Xn, yi, 2/2; • ■ ■ ) Un}), so that the result returned by this invocation to 
the algorithm is the best possible recombination. 



Random Column 
i 



Random 
Row 





A2 


As 


A4 



Bi 


B2 




Ai 


B2 


B3 


Ba 




B3 


A4 



Figure 6: Blind recombination operator for the MDSLP. 

In order to analyze the time complexity for this recombination operator, the critical part 
of the algorithm is the execution of lines 4-8 in Fig. 5. In this case, line 6 has complexity 
O(n^) (finding the minimum of at most 2n alternatives, the computation of each being 
0(n)). Line 6 has to be executed n/2 x 2n x 2n times at most (recall the optimization), 
making a global complexity of 0{7i^) = 0(|xp'^), where G 0(n^) is the size of solutions. 
Notice also that the recombination procedure can be readily made to further exploit the 
symmetry of the problem, extending variable domains to column values in addition to row 
values. The complexity bounds remain the same in this case. 

One interesting property of the described operator is that it can be generalized to re- 
combine any number of board configurations like BE-Opt(n, Uxesi-*"* M € {1 • • ^}})) where 
S" is a set comprising the solutions to be recombined. In this situation, the time complexity 
is 0{k^n^) (line 6 is 0{k'n?), and it is executed 0{k'^rfi) times), where k = l^l is the number 
of configurations being recombined. This multi-parental capability will be explored in the 
rest of the paper. 
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4.1.4 Experimental Results 

In order to evaluate the usefulness of the described hybrid recombination operator, a set 
of experiments for problem sizes from n = 12 up to n = 20 has been realized (recall that 
optimal solutions to the MDSLP are known up to n = 20). The experiments were performed 
using a steady-state evolutionary algorithm {popsize = 100, pm = l/ra^, px = 0.9, binary 
tournament selection). With the aim of maintaining diversity, duplicated individuals were 
not allowed in the population. Algorithms were run until an optimal solution was found 
or a time limit was exceeded. This time limit was set to 3 minutes for problem instances 
of size 12 and was gradually incremented by 60 seconds for each size increment. For each 
algorithm and each instance size, 20 independent executions were run. All the experiments 
in this paper have been performed in a Pentium IV PC (2400MHz and 512MB of main 
memory) under SuSE Linux. 

The base algorithm used is a MA using 2D-SPX for recombination, and endowed with 
tabu search for local improvement (maxiter = 11?). This algorithm is termed MAts, and 
has been shown to be capable of finding feasible solutions systematically, solving to op- 
timality instances with n < 15 (see MAts in Fig- 7). Although the performance of the 
algorithm degrades for larger instances, it provides distributions for the solutions whose 
average relative distance to the optimum is less than 5.29% in all cases. This contrasts with 
the case of plain EAs, which are incapable of finding even a feasible solution in most runs 
(Gallardo et al., 2006a). 



10 




instance size 



Figure 7: Relative distances to optimum for different algorithms for sizes ranging from 12 up 
to 20. Each box summarizes 20 runs. In this and in all subsequent figures, boxes comprise 
the second and third quartiles of the distribution (i.e., the inner 50%), an horizontal line 
marks the median, a plus sign indicates the mean, and circles indicate results further from 
the median than 1.5 times the interquartile-distance. 

MAts is firstly compared with MAs endowed with BE for performing recombination. 
Since the use of BE for recombination has a higher computational cost than a simple blind 
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recombination, and there is no guarantee that recombining two infeasible solutions will 
result in a feasible solution, we have defined three variants of the MAs: 

• In the first one, denoted MA-BE, BE is always used to perform recombination. 

• In the second one, termed MA-BEip, we require that at least one of the parents is 
feasible in order to apply BE; otherwise blind recombination is used. 

• In the last one, identified as MA-BE2F, we require the two parents to be feasible, thus 
being more restrictive in the application of BE. 

By evaluating these variants, we intend to explore the computational tradeoffs involved in 
the application of BE as an embedded component of the MA. For these algorithms, mutation 
was performed prior to recombination in order to take advantage of good solutions provided 
by BE. Fig. 7 shows the empirical performance of the different algorithms evaluated (as 
the relative distance to the optimum). Results show that MA-BE improves significantly 
over MAts and can find better solutions. MA-BE2F can find slightly better solutions 
than MA-BE on smaller instances (n G {13, 15, 16}), but on larger instances the winner is 
MA-BE. It seems that the effort saved not recombining unfeasible solutions does not further 
improve the performance of the algorithm. Note also that, for larger instances, MA-BEif is 
better than MA-BE2F- This correlates well with the fact that BE is used more frequently 
in the former than in the latter. 

As mentioned in Section 4.1.3, the optimal recombination scheme we use can be readily 
extended to multi-parent recombination (Eiben, Raue, & Ruttkay, 1994): an arbitrary 
number of solutions can contribute their constituent rows for constructing a new solution. 
Additional experiments were done to explore the effect of this capability of MA-BE. Fig. 8 
shows the results obtained by MA-BE for a different number of parents being recombined 
(arities 2, 4, 8 and 16). For arity = 2, the algorithm was able to find the optimum solution 
for all instances except for n = 18 and n = 20 (the relative distance to the optimum for the 
best solution found is less than 1.04% in these cases). Executions with arity = 4 cannot 
find optimum solutions for the remaining instances, but note that the distribution improves 
in some cases. Clearly, the performance of the algorithm deteriorates when combining more 
than 4 parents due to the higher computational cost. Variable clustering could be used to 
alleviate this higher computational cost, but this results in performance degradation since 
the more coarse granularity of the information pieces hinders information mixing (Cotta & 
Troya, 2000; Gahardo et al., 2006b). 

4.2 A Beam Search/MA Hybrid Algorithm 

In this subsection, we describe a hybrid tree search/memetic algorithm for the MDSLP. This 
algorithm combines, in a collaborative way, a BS algorithm and a MA. As noted before, 
BS works by extending in parallel a set of different partial solutions in several possible 
ways, and thus can be used to provide promising partial solutions to a population based 
search method such as a MA. The goal is to exploit the capability of BS for identifying 
probably good regions of the search space, and the strength of the MA for exploring these, 
synergistically combining these two different approaches. 

The proposed hybrid algorithm, that executes BS and the MA in an interleaved way, is 
depicted in Fig. 9. In the pseudo-code, a (possible partial) solution for an n x n instance 
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Arity=2 




instance size 



Figure 8: Relative distances to optimum for different arities for MA-BE for sizes ranging 
from 12 up to 20. Each box summarizes 20 runs. 

is represented by a vector of rows s = (ri,r2, . . . ,rj), i ^ n, where rows are encoded as 
binary strings, s • (rj = v) stands for the extension of partial solution s by assigning value 
V to its i-th row, and v denotes the reflection value of the binary string v. The hybrid 
algorithm, Hybrid(n, /c;,^, /cma), constructs a search tree, such that its leaves consist of all 
possible board configurations of size n x n that can be generated using solely symmetric 
rows (this symmetry constraint was imposed to keep the branching factor, kext, of the BS 
at a manageable level for the range of instance sizes considered), and internal nodes at 
level i represent partially specified (up to the i-th row) board configurations. This tree 
is incompletely traversed in a breadth first way using a BS algorithm with beam width 
fcft^ (i.e., maintaining only the best ki,^ nodes at each level of the tree). For the beam 
selection (line 10), a simple quality measure is defined for partial solutions, whose value is 
either oo if the partial configuration is unstable, or its number of dead cells otherwise. The 
algorithm starts (line 2) with a totally unspecified solution (i.e., a solution with rows). 
Initially, only the BS part of the algorithm is executed. During each iteration of the BS 
(lines 3-17), a new row is added to every solution in the beam (line 7). The interleaved 
execution of the MA starts only when partial solutions in the beam have at least kMA rows 
(line 11). For each iteration of the BS, the best popsize solutions in the beam are selected 
(using the quality measure described above) to initialize the population of the MA (line 12). 
Since these are partial solutions, they must be first converted into full solutions, e.g., by 
completing remaining rows randomly. After running the MA, its solution is used to update 
the incumbent solution (sol), and this process is repeated until the search tree is exhausted. 

4.2.1 Experimental Results 

Experiments were conducted to evaluate the hybrid algorithm (BS-MA-BE). The method- 
ology was the same as Section 4.1.4 (20 executions were performed for each algorithm 
and instance size), but arities for the MA where in {2,3,4}. The setting of parame- 
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Hybrid algorithm for the MDSLP 



function Hybrid {n, kh^, kMA) 

1 : sol := oo 

2: B :={{)} 

3 : for z := 1 to n do 

4: B':={} 

5 : for s ^ B do 

6: for r := to 2^^/21 - 1 do 

7: B' := B' U {s ■ in = r or r)} 

8: end for 

9: end for 

10: B := select best kbw nodes from B' 

11: if (i > k^fji) then 

12 : initialize MA population with best popsize nodes from B' 

13 : run MA 

14: sol := min (so/, MA solution) 

15 : end if 

16 : end for 

17: return sol 

end function 



Figure 9: Hybrid algorithm for the MDSLP. 



ters was kbw = 2000 (preliminary tests indicated that this value was reasonable), and 
kMA G {0.3 • n, 0.5 • n, 0.75 • n}, i.e., the best 2000 nodes were kept on each level of the BS 
algorithm, and 30%, 50% or 75% of the levels of the BS tree were initially descended before 
starting to run the MA. With respect to termination conditions, each execution of the MA 
within the hybrid algorithm consisted of 1000 generations, and no time limits were imposed 
for the hybrid algorithms that were run for n iterations of the BS. 

Fig. 10 shows the results for different values of parameter kMA- In order to better 
compare the distributions, the number of optimal solutions obtained by each algorithm 
(out of 20 executions) is shown above each box plot. For kMA = 0.3 • n, the performance of 
the resulting algorithm improves significantly over the original MA. Note that BS-MA-BE, 
using an arity of 2 parents, is able to find the optimum for all cases except for n = 18 
(this instance is solved with arity = 4). All distributions for diff'erent instance sizes are 
significantly improved. For n < 17 and arity S {2,3,4}, the algorithm consistently finds 
the optimum in all runs. For other instances, the solution provided by the algorithm is 
always within a 1.05% of the optimum, except for n = 18, for which the relative distance 
to the optimum for the worst solution is 1.3%. The other two charts show that, in general, 
the performance of the algorithm deteriorates with increasing values of the kMA parameter. 
This may be due to the low quality of the bounds used in the BS part. 
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Figure 10: Relative distances to optimum for different arities for BS-MA-BE and Kma £ 
{0.3 ■ n,0.5 • n,0.75 • n}, for sizes ranging from 12 up to 20. Each box summarizes 20 runs. 
The numbers above each box indicate how many times the optimal solution was found. 



Regarding execution times, Fig. 11 shows the distributions for the time (in seconds) 
to reach the best solution needed by the algorithms. Although BS-MA-BE requires more 
time than MA-BE, the time needed remains reasonable for these instances, and is always 
less than 2000 seconds. Note also how the execution time increases with the arity, as more 
time is needed by the MA to perform BE in the crossover operator. On the other hand, 
execution time decreases for larger values of kMA as the number of executions of the MA 
decreases, although, as we have already remarked, the quality of the solutions worsens. 



To verify that the improved results of the hybrid algorithm were not only a consequence 
of the extended execution times, experiments for MA-BE were repeated with an increased 
time limit of 2800 seconds for each execution independently of the instance size. The results 
of these experiments are shown in Fig. 12. Clearly, the performance of MA-BE does not 
improve dramatically, and this provides evidence on the synergetic cooperation of BS and 
MA achieved by the hybrid algorithm. 
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Figure 11: Time (in seconds) to best solution for different arities for BS-MA-BE and Km A £ 
{0.3 • n, 0.5 • n, 0.75 • n}, for sizes ranging from 12 up to 20. Each box summarizes 20 runs. 



5. A New Hybrid Algorithm Based on Mini-Buckets 

In this section we present a novel hybrid algorithm based on the algorithm described in 
Section 4.2. This algorithm exploits the technique of Mini-buckets that is explained in the 
following. 

5.1 Mini-Buckets 

The main drawback of BE is that it requires exponential space to store functions extension- 
ally. When this complexity is too high, the solution can be approximated using the technique 
of mini-buckets (MB) presented by Dechter (1997) (see also Detcher & Rish, 2003). Recall 
that, in order to eliminate variable Xj, with its corresponding bucket Bi = {fi^, ■ ■ ■ fi^}, BE 
calculates a new cost function 

5i = ( ^ /) ^ (27) 
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Figure 12: Relative distances to optimum for different arities for MA-BE executed for 2800 
seconds, for sizes ranging from 12 up to 20. Each box summarizes 20 runs. 



whose time and space complexity increases with the arity of gi, i.e., with the arity of the 
set U/GSi ^"^^(Z) ~ {xi}- This complexity can be decreased by approximating the function 
Qi with a set of smaller- arity functions. The basic idea is to partition bucket Bi into k so 
called mini-buckets Bi-^ , ■ ■ ■ , Bii^, such that the number of variables in the scope of each Bi. 
is bounded by a parameter. Afterwards, a set of k cost functions with the reduced arity 
sought can be defined as 

9i, = iYl f)^x^,j = 'i^---k, (28) 
and the required approximation to gi can be computed as their sum: 

k k 

ft = E^^^. = E ((E (29) 

Note that the minimization computed in gi by the ij. operator has been migrated inside the 
sum. Since, in general, for any two non- negative functions fi{x) and f2{x), minx{fi{x) + 
f2{x)) > minxfi{x) +minxf2{x), the following inequality holds 

gi r-^ V 

(E/)^^''^E((E (30) 

/GSi j=l /GBi. 

and, thus g\ is a lower bound on gi. Therefore, if variable elimination is performed us- 
ing approximated cost functions, it provides a lower bound for the optimal cost requir- 
ing less computation than BE. Notice that the described approach provides a family of 
under-estimating heuristic functions whose complexity and accuracy is parameterized by 
the maximum number of variables allowed in each mini-bucket. 
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5.2 Improving the Lower Bound Using Mini-Buckets 

The simple quahty measure for beam selection used in the algorithm in section 4.2 depends 
solely on the part of the solution that is already constructed. In this section, we will use 
the MB technique to compute a tight, yet computationally inexpensive, lower bound for 
the remanning part of the configuration with the aim of improving the performance of the 
BS part of the hybrid algorithm. 

For this purpose, let us note that the MDSLP for an n x n board can be formulated 
as an alternative WCSP, if we associate a different variable Xij for each cell on the 
board. With this formulation, there are cost functions fij,l ^ i,j ^ n. The scope of 
function fij is Xij and all its neighborhood, and it returns oo if the cell is unstable, 1 
if cell is dead, and otherwise. The following objective function 

n n 

^ = EE/«^- (31) 
i=i j=i 

has to be minimized. 

Note that the original formulation, introduced in Section 3.2, can be obtained from the 
present one by clustering all cost functions corresponding to row rj into a single cost function 
fi. In the same way, let us cluster cost functions for the i-th row into M cost functions, 
fl, ff, . . . , f^^ of roughly the same arity n/M), each one evaluating respectively one of 
the M segments of the row. To be precise, 

Em 

/r= E 1^"^^^^ (32) 

where Wn, 1 ^ n ^ M, stands for the number of variables of each segment. 

Using this formulation, BE would perform the elimination of all variables corresponding 
to the last row by computing a new gn cost function as 

M M 

9n = {Y.fn-l+Y.fn)^{^nl:Xn2,---Xnn}, (33) 
m=l m=l 

whose bucket is Bn = • • • , /«. /«'•••. fn}- Applying mini-buckets, S„ 

can be partitioned into M buckets: -B™ = < m < M, and a set of M cost 

functions to approximate gn with reduced arity can be calculated as: 

gn = ifn-l + fn)i^xT, l<m<M (34) 

where 

(35) 

'^i(i+E}1i '^<m<M (36) 

iXin\- (37) 



M _ r 
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In this way, the number of variables in each meta-variable x^, 1 < m < M, is n/M 
approximately. Because the scopes of g^, 1 < m < M, are {x^_2,x^_^}, their arities are 
approximately 1/M of the arity of gn- The rest of the rows of the board can be processed 
in a similar way. 

Cost functions computed by the function MB can be used to estimate a tight lower 
bound for a partial solution during the execution of the hybrid algorithm as follows: let 
s = (ri, r2, . . . , Vk) be a partial solution with k rows for an n x n instance of the MDSLP. 
As defined, gl^i{rl_^,rl) returns the cost of the best extension to partial solution s that 
can be attained in rows k to n, considering only the first column. In a similar manner, 
,1 < m < n can be used to estimate the best extension considering only columns 2 to 
n respectively. Hence, a lower bound for a partial solution can be computed as: 

fc— 1 n n 

lbin,r,, . . . , rfc) = ^ + ^ <7r+i(^r-i, ^D, (38) 

i=l j=l m=l 

where the first sum corresponds to the part of the solution already assigned. This bound 
can be used to rank nodes for beam selection and the initialization of the MA population. 
In the following subsection, we have experimented with setting M = 3, so that 



' (n/3, n/3, n/3), n mod 3 = 

([n/3j, [n/3], [n/3j), n mod 3 = 1 (39) 
J[n/3], [n/3j, [n/3]), n mod 3 = 2. 



Observe that, for these settings, the space complexity of function MB is 0{n x 2'^^^a~\~^'^^), 
whereas its time complexity is 0(n^ x 2^''I^3"T+^)). When this complexity is still too high, 
the approach described in this subsection can be utilized to reduce it further, considering 
more than three clustered cost functions for each row of variables, although the resulting 
bounds would be less tight. 



5.2.1 Experimental Results 

Experiments were repeated for the hybrid algorithm equipped with the new lower bound, 
BS-MA-BE-MB. Fig. 13 shows the results of these experiments for values of kjy.jA £ {0.3 • 
n, 0.5 • n, 0.75 • n}. The algorithm finds the optimum for all instances and arities and the 
relative distance to the optimum for the worst solution found is less than 1.05% in all cases. 
The best results are obtained with arity = 4, although this requires slightly more execution 
time. Note also how BS-MA-BE-MB is less sensitive to the setting of parameter kMA, 
which means that execution times can be reduced considerably using a large value for this 
parameter (see Fig. 14). The particular combination of parameters kMA = 0.75 • n and 
arity = 4 provides excellent results at a lower computational cost, as execution times are 
always below 570 seconds for n ^ 20. As a comparison, recall that the only approach in the 
literature that can solve these instances - described by Larrosa et al. (2005) - requires over 
33 minutes for n = 18, 15 hours for n = 19 and 2 days for n = 20, and that other approaches 
are unaffordable for n > 15. Note however that these times correspond to a computational 
platform different to ours. In order to do a fairer comparison, we executed the algorithm of 
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Figure 13: Relative distances to optimum for different arities for BS-MA-BE-MB and 
Kma £ {0.3 • n,0.5 • n,0.75 • n}, for sizes ranging from 12 up to 20. Each box summa- 
rizes 20 runs. 
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Figure 14: Time (in seconds) to best solution for different arities for BS-MA-BE-MB and 
kMA = 0.75 • n, for sizes ranging from 12 up to 20. Each box summarizes 20 runs. 



Larrosa et al. ^ in our platform. In this case, it required 1867 seconds (i.e., more than 31 
minutes) in order to solve the n = 18 instance, and more than 1 day and 18 hours to solve 
the n = 20 instance. These values are very close to the times reported by Larrosa et al. 
(2005), and hence indicate that the computational platforms are fairly comparable. 

5. Available at http : //www. Isi .upc . edu/~larrosa/publications/LIFE-SOURCE-CODE .tar . gz . Time for 
n = 19 could not be obtained as the code provided by Larrosa et al. can only be used with even sized 
instances. 
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□ Arity=4 
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Figure 15: Relative distances to best known solutions for different arities for 
BS-MA-BE-MB and /cam = 0.3 -n, for very large instances (i.e., sizes of 22, 24, 26, and 28). 
Each box summarizes 20 runs. Note the improvement of best known solutions for sizes 24 
and 26. 
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Figure 16: New best known maximum density still lifes for n € {24,26}. 



Table 2: Optimal solutions for the SMDLP. 
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opt 


232 
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378 













5.3 Results on Very Large Instances 

As already mentioned, there is currently no approach available to tackle the MDSLP for 
n > 20. Larrosa et al. (2005) tried their algorithm for n = 21 and n = 22, but they could 
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not solve any of those instances within a week of CPU). For these very large instances, only 
solutions to some relaxations of the problem are known. One of these relaxations, known 
as the symmetrical maximum density still life problem (SMDSLP), was proposed in (Bosch 
&; Trick, 2002), and consists of considering only symmetric boards (either horizontally or 
vertically) which reduces the search space from 2" to 2" ^'^1'^ . 

The optimized version of BE algorithm (Section 2.4) can be used find vertically sym- 
metric still lifes, by defining as the variable domain, a set that contains only symmetric 
values for rows, 

p = {r or r I r € {0..2r"/2l _ 1}}. (40) 

Larrosa and Morancho (2003) and Larrosa et al. (2005) used this algorithm to solve the 
SMDSLP for the instances considered so far in this paper (i.e., for n € {12 . . 20}), as well as 
for very large instances (i.e., n € {22, 24, 26, 28}). Results are summarized in Table 2, which 
shows for each instance size the optimal symmetrical solution (as the number of dead cells). 
Clearly, the costs of optimal symmetric still lifes are upper bounds for the MDSLP, that 
can additionally be observed to be very tight for n ^ 20. Results for n > 20 are currently 
the best known solutions for these instances. 

We also run our algorithm (BS-MA-BE-MB) for these very large instances (i.e., n € 
{22,24,26,28}), and compare our results to symmetrical solutions for these instances. Re- 
sults (displayed in Fig. 15 shows that our algorithm was able to find two new best known 
solutions for the MDSLP, namely for n = 24 and n = 26. There are 275 and 324 dead 
cells respectively in the new solutions. These solutions are pictured in Fig. 16. Incidentally, 
our algorithm could also find a solution with 325 dead cells for the n = 26 instance. For 
the other instances, our algorithm could reach the best known solutions consistently. Let 
us note that the computation of mini-Buckets for these very large instances was done by 
considering four clustered costs functions for variables in each row of the board, as the 
complexity when using three costs functions was still too high. 

6. Conclusions 

The MDSLP represents an excellent example of WCSP; its highly constrained nature is typ- 
ical in many optimization scenarios. Furthermore, the algorithmic hardness of solving this 
problem illustrates the limitations of classical optimization approaches. For this reason, it 
is not surprising that this problem has attracted the interest of the constraint-programming 
community, and has been central in the development and assessment of sophisticated tech- 
niques such as bucket elimination (BE). However, the high space complexity of BE as an 
exact technique (Dechter, 1999), makes this approach impractical for large instances. In 
this work, we have presented several proposals for the hybridization of BE with memetic 
algorithms and beam search (BS), and showed that they represent very promising models. 
The experimental results have been very positive, solving to optimality large instances of 
the MDSLP. We have also studied the influence that multi-parent recombination have on 
the performance of the algorithm. The results indicate multi-parent recombination can help 
to improve the results obtained by previous approaches. 

Among all our proposals, we must distinguish a new algorithm resulting from the hy- 
bridization, at different levels, of complete solving techniques (i.e., bucket elimination), 
incomplete deterministic methods (i.e., beam search and mini-buckets) and stochastic algo- 
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rithms (i.e., memetic algorithms), that empiricahy produces good-quahty results, not only 
solving to optimality very large instances of the constrained problem in a relatively short 
time, but also providing new best known solutions in some large instances. This algorithm 
exploits the technique of the mini-buckets to compute tight yet computationally inexpen- 
sive lower bounds of the partial solutions that are considered in the BS part of the hybrid 
algorithm. 

As future work, we plan to consider complete versions of the hybrid algorithm. This 
involves the use of adequate data structures to store not yet considered but promising 
branch-and-bound nodes. While the memory requirements will of course grow enormously 
with the size of the problem instance considered, it will be interesting to analyze the com- 
putational tradeoffs of the algorithm as an anytime technique. 
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